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Abstract 

Three problems in the statistical mechanics of models for an assembly of molecular 
motors interacting with cytoskeletal filaments are reviewed. First, a description of 
the hydro dynamical behaviour of density-density correlations in fluctuating ratchet 
models for interacting molecular motors is outlined. Numerical evidence indicates 
that the scaling properties of dynamical behavior in such models belong to the KPZ 
universality class. Second, the generalization of such models to include boundary 
injection and removal of motors is provided. In common with known results for the 
asymmetric exclusion processes, simulations indicate that such models exhibit sharp 
boundary driven phase transitions in the thermodynamic limit. In the third part of 
this paper, recent progress towards a continuum description of pattern formation in 
mixtures of motors and microtubules is described, and a non-equilibrium "phase- 
diagram" for such systems discussed. 



1 Introduction 



Living systems exhibit a remarkable variety of non-equilibrium steady states. 
Problems associated with the modeUing of such states include the description 
of the non-equilibrium behaviour of membranes driven by active pumps[l], hy- 
drodynamic approaches to the motion of self-propelled objects[2,3], the theory 
of pattern formation in a variety of biological contexts[4] and models for intra- 
cellular transport processes associated with the motion of molecular motors[5]. 
This last problem has attracted the attention of statistical physicists in recent 
years, since the simplest models for such systems have several advantages: they 
are exactly solvable even in the presence of interactions, easy to generalize. 
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relatively straightforward to simulate and closely related to models studied 
extensively in the context of traffic fiow[5]. 

It is useful to develop intuition with simple models, requiring only that it 
should be possible to add the requisite biological detail incrementally. This 
would then allow progressively more accurate descriptions to be constructed 
within a sequence of increasingly refined models. This paper reviews some cal- 
culations which explore the middle ground between extremely simplified statis- 
tical mechanics models for the motion of interacting molecular motor proteins 
- the asymmetric exclusion process and generalizations - and marginally more 
realistic models for the motion of individual motor proteins generalized to ac- 
comodate motor-motor interactions [6]. It also reviews some recent work on 
the hydrodynamic description of pattern formation in mixtures of molecular 
motors and microtubules[7]. 

Molecular motors are a class of biological machines which function within 
cells[8]. Such motors, proteins such as kinesins, myosins and dyneins, move 
unidirectionally on one-dimensional "tracks" while hydrolysing adenosine tri- 
phosphate (ATP). These tracks are components of the cytoskeleton in eukary- 
otic cells, an extended dynamic network formed through the polymerization 
and crosslinking of tubulin and actin monomers to form microtubules and 
actin filaments[9]. The cytoskeleton helps the cell anchor to substrates, to 
move and to divide, and lends it mechanical and structural rigidity. In addi- 
tion, this network defines paths for molecular motors to transport cargo to 
different parts of the cell [9]. 

The asymmetric nature of motor motion along the cytoskeleton derives from 
the asymmetry of the constituent monomeric units of the track. Microtubules 
( equivalent ly, actin filaments) can be idealised as periodic, one-dimensional, 
rigid structures, this periodicity following from their polymeric nature. The 
violation of detailed balance necessary in order for the motor to exhibit di- 
rected motion comes from the transduction of the chemical energy obtained 
from ATP hydrolysis into mechanical work[10]. The stochastic uptake of ATP 
is one source of random noise in the problem; the other is the thermal noise 
which dominates all biological systems at cellular scales. Molecular motors 
thus act as Brownian rectifiers in exhibiting a non-zero drift velocity in the 
absence of a net time-averaged force. The problem of modelling molecular mo- 
tors can therefore be placed in the more general context of "Brownian motor" 
or "thermal ratchet" models for the extraction of useful work from thermal 
fiuctuations[10,ll]. 

We consider the "fiuctuating potential" or "fiashing ratchet" model for Brow- 
nian motors[10]. In this model, a single motor, idealized as a point object 
moving in one dimension, is driven by stochastic forces uncorrelated in space 
and time and drawn from a Gaussian distribution with zero mean. The motor 
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Fig. 1. The continuum potential corresponding to the two states (non-trivial and 
flat) of the potential, showing the locations of several motors interacting through a 
hard-core repulsive potential. 

also experiences a force derived from the gradient of a time-dependent poten- 
tial. The time dependence of this potential is generated by switching randomly 
between two states, one in which the potential is an asymmetric sawtooth - 
say with the longer leg of the potential to the left - and the other in which it 
is flat. 

In the "off' or "fiat" state of the potential, particles diffuse isotropically. When 
the potential is switched to the "on" state, a particle is more likely to be found 
in a region of space where it experiences a net force to the left than to the right, 
given that the potential lacks right-left symmetry. Repeated cycling between 
"off" and "on" states generates net motion. The generation of directed motion 
occurs through a subtle mechanism: in any one of the potential states, given 
sufficient time to equilibrate, no net current can flow provided the microscopic 
jump rates obey detailed balance. The breaking of detailed balance overall 
arises from the non-equilibrium, time-dependent switching between potential 
states and not from the choice of hopping rates in any one of these states [10]. 

In modelling biological motors, the asymmetric potential encodes the energy 
associated with an internal state of the motor as a function of its position 
along the track[10]. This state represents a particular conformation of the mo- 
tor protein. ATP consumption induces transitions between states. A motor 
molecule thus has two degrees of freedom, a spatial coordinate and an inter- 
nal (state) coordinate. A faithful representation of the internal states of the 
motor and track thus involves specifying a large number of continuous peri- 
odic potentials representing the energy of a molecular motor at location x as 
a function of its state s. 

The simplified model described here retains only two states of the motor- 
track complex. Such "minimal" models may be expected to capture some of 
the relevant complexity of the real biological system, at least in the limited 



contexts in which we will study them. It is straightforward to generalize the 
models described here. In particular, one could allow motors both to detach 
from the filament and diffuse within the ambient solvent as well as to reattach 
with prescribed rates. One could also allow for motors of finite extent, by 
extending the exclusion constraint to sites which neighbour the one occupied 
by the motor. 

The effects of cooperativity between motors has been a source of rich physics 
in recent years [10,12,13,14,15]. Such effects, in biological models for the action 
of the myosin motors involved in muscle contraction, have been incorporated 
by coupling motors through intervening elastic elements. Thus, in addition to 
forces derived from the external potential and Gaussian noise, a motor feels a 
force due to its elastic interactions with its neighbours. A mean-field analysis 
of the effects of such elastic couplings yields a variety of novel phenomena, such 
as Hopf bifurcations and spontaneous oscillations of the current [10]. However, 
another class of interaction effects which operate between elastically decoupled 
motors can be envisaged: the steric hindrance (or, in general, any short-range 
interaction) experienced by motors translocating on a filament when they 
approach each other [16, 6, 17]. 

How can the fluctuating potential model be generalized to incorporate the 
effects of interactions between many motors moving on a cytoskeletal filament? 
The simplest such interaction is a hard-core interaction between two motors, 
as illustrated schematically in Fig 1. On the lattice, this constraint is simply 
implemented by allowing only one particle to occupy a given lattice site at 
a time. We define simple lattice models incorporating such interactions and 
study these models with boundary conditions (periodic) which conserve the 
total number of motors in Section II. We also discuss a simple mean-field 
theory and its prediction for currents and density profiles within a single period 
of the potential. A second set of results, presented in Section HI, relate to 
the existence of phase transitions in such systems induced by the effects of 
adding (subtracting) motors at the boundaries of the (open) one-dimensional 
chain at prescribed rates [18]. Our numerically calculated phase diagram closely 
resembles the phase diagram of the partially asymmetric exclusion process 
with boundary driving[19,20]. 

The fact that these two systems should be related at the level of hydrody- 
namic correlations was conjectured in Ref.[6] and used to predict that density- 
density correlations in the steady state of the interacting motor system should 
obey scaling with exponents belonging to the KPZ universality class. The 
fact that these systems show the same type of phase transitions as a func- 
tion of boundary conditions is further evidence of the close relationship be- 
tween these models, despite the far greater complexity (multiplicative time 
and space-dependent noise at the level of the microscopic hopping rates) in 
the Brownian motor model. 



A third set of results reviewed in this paper (Section IV) relate to the mod- 
elling of the patterns which form when motor complexes are mixed with mi- 
crotubules and supplied with ATP, in a quasi- two-dimensional geometry [7]. 
These patterns include structures such as asters, vortices and aster-vortex 
mixtures as well as disordered states. Understanding the generic features of 
these states, the sequence of transformations between them as a function of 
the motor density and the interactions which contribute to the formation of 
such self-organized structures, is believed to be a crucial part of understand- 
ing the physics behind the formation of a cellular-scale pattern universal to 
all eukaryotic cells, the mitotic spindle [21,22,23,24,25,26,27]. The concluding 
section. Section V, discusses some general features of the results, suggesting 
that the general attribute of "physical robustness" , a robustness of the non- 
equilibrium steady states obtained in such models towards a wide class of 
physically relevant perturbations, may be biologically relevant. 



2 KPZ Scaling in Interacting Ratchet Models 



In the simplest continuum versions of the fluctuating potential models, indi- 
vidual motor particles see a time-dependent potential V{x,t) = r](t)U{x), in 
addition to random Brownian forces with zero mean value. Here U{x) is peri- 
odic with period £ i.e. U{x + i) = U{x) and an asymmetric function of x i.e. 
U{x) 7^ U{—x). The time dependence of V{x,t) is governed by a (stochastic 
or deterministic) switching function 77 (t) which takes the values and 1. We 
assume U{x) to be of the sawtooth form 

U{x)=ax {0<x<uji), 

= b{i-x) {iui<x<i) (1) 

with a, b > 0, aui = bi{l — uj) and uj < 1. The switching of the potential 
occurs independently of the state of the motors, thus breaking detailed bal- 
ance. Together with the lack of reflection symmetry in U{x), this switching 
generates a net particle current. 

A useful simplification is the discretization of the periodic potential in space to 
convert the continuum problem into a lattice one. (This is not a unreasonable 
simplification, since the biological motor appears to undergo a sequence of 
discrete conformational changes, each of which is coupled to a partial translo- 
cation across the period.) Each period of the potential is divided into W 
lattice sites, all of which are assigned to the segment of the potential with 
positive slope. The length of the system, L, is measured in periods of the 
sawtooth. The maximum height of the sawtooth potential is V and we de- 
fine r = exp(— ^ TOy-i) )- Filially, the parameters Pqi and Piq represent the 
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Fig. 2. Dependence of the steady-state current on the density p of particles for 
32 < L < 256. 

probabilities that the potential changes from 77 = to 77 = 1 and vice versa. 

The transition rates between the configurations of motors are chosen to satisfy 
detailed balance with Metropohs rates, P ({c'} -^ {o"}) = min(l, exp[iJ({(T'}) — 
H[{a})\). Here {a} indexes allowed configurations of the motors on the lattice 
and H{{(j)} is the energy of the configuration. In terms of these parameters 
the hopping probabilities for motor particles on the sawtooth are 
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where we have indexed the lattice sites from at the potential minimum. The 
hopping probabilities for left(right) jumps on the flat potential equal 1/2 at 
every site. 
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Fig. 3. The relaxation function $(A;, t) for L = 128 and p = 0.5 plotted as a function 
of the scaled variable k^t for z = 1.60 for the five smallest values of k = lirj/L, 
with j = 1,2,3 .. .. The lowest k- value splits off whereas data for all higher k fall on 
the same branch. 



To generalize models for single Brownian motors to a finite number Np of 
interacting motors on a one-dimensional lattice, we assume that these mo- 
tors hop to unoccupied nearest neighbour sites with rates determined by the 
discretized potential in Eq. (2), exactly as they would in the non-interacting 
case. The only interaction between these motors is thus a hard-core repulsion 
which prevents them from occupying the same lattice site. Our discretization 
and the hard-core constraint ensures an upper bound on the number of motors 
which can occupy the lattice. 



We consider periodic boundary conditions in the calculations described in 
this section. Our motors thus move on a ring. An elementary step consists of 
either an attempt of a particle to hop to a neighbouring site or an attempt to 
switch the potential state globally. The results of Ref. [6] were obtained by a 
procedure which involved changing the potential state globally. Flipping this 
state locally does not alter the conclusions qualitatively and the quantitative 
changes are small. 
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Fig. 4. The effective exponent Zeff{N, p) plotted as function of p for systems of size 
32 < L < 256 Note the systematic decrease of Ze// ^is L is increased at constant 
density. 

Numerical results for the model defined above (Model I of Ref. [6]) are obtained 
in the following way: We use W = 10 lattice sites per period of the asymmetric 
sawtooth in all our simulations, varying the system size L = 24,48,64,128, 256 
and 512. We took r = 0.05, Pqi = 0.03 and Pw = 0.04, where Pqi and Pw 
are the probabilities that the potential state goes from flat to non-flat and 
vice versa. Typically, the system equilibrates over 5 x 10^ MCS before data for 
currents and correlation functions are recorded. These quantities are averaged 
over 5 x 10^ — 10^ conflgurations. 

Fig. 2 illustrates the fundamental relation (the functional dependence of the 
current on the density) of the model. We examine principally the scaling prop- 
erties of the intermediate scattering function Sp{k,t) deflned by 



S,{k,t) = ^<Sp{k,0)Sp{-k,t)>, 



(3) 



where A^ is the number of particles, k is 2mT/L, n = 0, ±1, ±2, . . . , L/2, and 
L is the system size. Here Sp{k, t) is the Fourier component with wave vector 



k of the deviation from the mean local density, Sp{x,t) = p{x,t)— < p{x) > 
with the brackets < • > denoting a time average. 



Transforming the density field p{x, t) to a "height" field h{x, t) via p{x, t) = 
dxh{x,t), and imposing helical boundary conditions on h{x,t) to satisfy the 
constraint /q p{x, t) = h{L) — h{0) = N, Sp{k, t) can be related to the structure 
factor S{k,t) = {6h{k,0)6h{—k,t)) where 6h{k,t) is the Fourier transform of 
h{x,t) — {h{x)). For small k and large t, i.e. in the hydrodynamic limit, if 
S{k,t) exhibits dynamical scaling, we can write S{k,t) ~ k~'^^^F{k^t), where 
1] and z are scaling exponents and F is a scaling function. 

In Ref.[6], it was conjectured that the dynamical scaling properties of the 
"height" field should be described by the Kardar-Parisi-Zhang [28] equation, 
a non-linear Langevin equation of the form, 

f = ^V^/^ + ^(V/.f + C(x,t), (4) 



which is known to describe the long-time, long- wavelength behavior of a num- 
ber of nonequilibrium systems [29]. This equation is written in a form appropri- 
ate for surface models where h{x,t) is the height, relative to a rf-dimensional 
substrate, of a growing interface and C(x;^) represents white noise. For the 
KPZ equation, owing to the existence of a fluctuation-dissipation theorem, 
the exponents can be obtained exactly for rf = 1 and take on the values 
z = 3/2 and rj = 0. The connection is made in the following way: Coarse-grain 
microscopic configurations of such models in space and time. At spatial scales 
larger than the repeat distance i of the periodic potential and for time scales 
much larger than the typical time-scale r over which the potential changes, the 
system will appear to have a constant density on average, as well as a constant 
current. Superimposed on this constant density are spontaneous fluctuations 
which obey a local conservation law. The effects of interparticle interactions 
at the largest length scales can be summarized in the following observation: 
These density fluctuations are convected with a speed, the "kinematic wave 
speed" , which depends on their magnitude. 

Consider now the statistical mechanics of an unrelated model, that of the 
stochastic dynamics of a flnite density p of hard-core particles on a line, which 
hop individually with rate (1 + e)/2 to the right and (1 — e)/2 to the left, pro- 
vided the excluded volume constraint is satisfled is known as the "asymmetric 
exclusion process" (ASEP) for e 7^ 0. The ASEP has a net particle current 
J = ep{l — p). The symmetry breaking which results in a constant current in 
the ASEP is an explicit consequence of the asymmetry in the hopping rates. 
This symmetry breaking is to be contrasted to the more subtle symmetry 
breaking in the case of the ratchet models. 



Density-density correlations for the ASEP are known to be governed by KPZ 
exponents[30,31,32]. Since both models share the feature expected to be most 
relevant to a hydrodynamic description - the existence of a non-trivial density- 
dependent current - it is reasonable to conjecture that they should belong to 
the same universality class, irrespective of the fact that the detailed origin of 
the symmetry breaking is different in each case. 

In our simulations we measured the relaxation function [30] 

^ <pi-k, 0)pik, t)>- <pi-k, 0) xpjk, t) > ^ Sp{k,t) 
^ ' ' ~ <p{-k, 0)p{k, 0) > - <p{-k, 0) ><p{k, 0) > Sp{k, 0) ' ^ ' 



which is just the Fourier transform of J2i{Pi(t) — Pi(0))(Pj+j(t) — Pi+j(0)), 
normalised by its value at t = 0, an arbitrarily chosen time in the steady 
state. The data for ^{k,t) for a given density and system size were plotted 
as a function of the scaled variable k^t for various z and the value of z that 
provided the best collapse of the data by visual inspection was taken to be Zeff- 
Examples of this data collapse are shown in Fig 3 at half filling for L = 128 
for Zeff = 1.60. The relaxation function generically has two distinct branches: 
For the smallest value oi k, k = 27r/L, the relaxation function decays more 
slowly than for larger values of k. For j > 1 the data collapse to quite high 
accuracy onto a single curve. This separation of the relaxation function into 
two branches (with j = 1 the special case) is also a feature of the single step 
model [33]. 

We have carried out this analysis systematically for L = 32 to L = 256 as 
function of p. The effective exponents Zeff{N, p) are plotted as functions of p in 
Fig. 4. There is a systematic decrease of the effective exponent Zeff as function 
of increasing L for all p ^ 0, with the smallest values occuring at p = 0.5, as 
expected on theoretical grounds. The smallest value of the critical exponent, 
z, from this data was found to be 1.58 ± .01 for L = 256. Extrapolating to the 
thermodynamic limit, z{p) approaches the function: 



^ for < p < 1 
2 for p = 0, 1. 

consistent with the predictions of the KPZ equation. 
2. 1 Mean-Field Approximation 



Assuming periodic boundary conditions, all periods of the potential are equiv- 
alent in steady state. By translational invariance, it is then sufficient to solve 
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for the steady state density fields and steady-state currents within a single pe- 
riod. If a site is to be updated, which is the case with probability 1/2 in each 
elementary time step, the density pi(t) at site i at time t changes through hops 
on and off that site from neighbouring sites. These hops can only occur if the 
hard-core constraint is satisfied. Allowed hops occur with the (microscopic) 
probabilities pf(t) and pf (t) for hops to the left (pf (t)) or right {pfit)) at site 
i at time t. Given our algorithm, in which an attempt is made to update either 
a particle or a potential state at each time step and never both, the assignment 
of particle hopping rates at time t is unambiguous. However, how these rates 
are to be interpreted in the time-continuum limit is not straightforward and 
will be dealt with explicitly in what follows. 

The update processes for a site i at time t can be written purely in terms of 
the local density fields in the neighbourhood of that site. The density field 
Pi{t + 1) at time t + 1 is given by : 



Probability (1 - 3/A^) 

Probability (1 -pf_-^)/N 
Pi{t))pi-i{t) Probability pf„i/N 

Probability (l- pf' - pf)/N 
Pi{t + l) = \ pi{t)pi_i{t) Probability pf/N 

Probability pf /N 

Probability (l-pf+i)/iV 

(6) 
piit))pi+i{t) Probability pf_^i/A^ 

The first term in Eq. (6) represents the probability that a site other than i is 
picked at time step t, whereas subsequent terms represent probabilities that 
the density variable at site i is updated as a consequence of site i — l,i or 
i + 1 being picked. An overall factor of 1/2 in these transition probabilities 
(arising from the fact that the choice to update a site or potential is made 
with probability 1/2 at every time step), has been absorbed into a rescaling 
of time. 

The hopping rates pf^_i, pf , pf , Pi^i are all exphcitly functions of time and 
are determined by the instantaneous state of the potential. This state is deter- 
mined by a corresponding equation for the dynamics of the potential field or 
equivalently of the stochastic variable r]{t) which appears in the definition of 
V{x,t). We must specify two averages in the steady state. One is the average 
over thermal noise given a particular stochastic history of potential fiips. The 
other is the average over all stochastic histories. We are interested in those 
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Fig. 5. The current-density relation as obtained through Monte Carlo simulations 
and the mean field theory described in the text, assuming r = 0.05, Pqi = 0.03 and 
Pio = 0.04. 

attributes of the system which characterize its steady state. The dynamics of 
potential and particle are partially decoupled. The potential state influences 
the particle hopping rate but the potential flips independently of the particle 
state. The mean-field approximation made is the following: Whenever a prod- 
uct such as pi(t)pi+i(t) is to be averaged over either thermal noise or potential 
flip histories (with the appropriate averaging operation denoted by (■)), we 
replace {pi(t)pi^i(t)) by (pi(t))(pj+i(t)). This approximation truncates the hi- 
erarchy of coupled correlation functions by representing correlation functions 
of all higher-order products of density fields in terms of single-site averages. 
Our mean- field theory is formulated for the following limit: If the potential 
fluctuates over a microscopic time-scale which is much faster than character- 
istic diffusion time scales over a single period, it is legitimate to average the 
rates. The discrete equation Eq. (6), can then be converted into a first-order 
non-linear system of differential equations in time. These equations are 



dpijt) 
dt 



rlpi{t) + rfpi.iit) + Tfpi+i{t) + Tfpi.i{t)pi{t) + T!pi{t)pi+i{t){7) 



where i runs over the sites in a single period. Given Pqi and Pio, 
of t} . . . rf are the following: r} 
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and Tj = Pi + Pi+i- The bars denote an average over the rates. The set of 
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Fig. 6. The comparison between local densities in each period as obtained by aver- 
aging the Monte Carlo data and predictions of the mean field theory described in 
the text, assuming r = 0.05 Pqi = 0.03 and Pio = 0.04. The mean density is 0.25. 



equations Eq.(7) represent the mean- field treatment of the case in which the 
average over stochastic histories of potential flips has been performed before 
the average over thermal histories. The steady state in mean-field theory is ob- 
tained by setting the time derivatives dpi{t)/dt to zero. The resulting equations 
are to be solved for the W sites within a period, Given a mean field solution 
for the densities, the time-averaged current in the mean-field approximation 
can be obtained. Fig 5 illustrates the comparison of the fundamental diagram 
of the system as obtained through the mean field theory with the Monte Carlo 
data. The density profiles for densities 0.25 and 0.75 are shown in Figs. 6 and 
7 together with the results of a direct numerical simulation. It can be seen 
that the mean-field theory in the averaged-rate limit used above predicts the 
currents and the density profiles to reasonable accuracy. At the level of the 
currents, the agreement between the mean-field theory and the simulation re- 
sults is certainly passable. The density profiles appear qualitatively accurate 
but are incorrectly rendered in quantitative terms. 
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Fig. 7. The comparison between local densities in each period as obtained by aver- 
aging the Monte Carlo data and predictions of the mean field theory described in 
the text, assuming r = 0.05 Pqi = 0.03 and Pio = 0.04. The mean density is 0.75. 

3 Boundary-Driven Phase Transitions 



The periodic boundary condition assumed in previous sections is clearly in- 
correct in the biological context, where motors are loaded and unloaded along 
the cytoskeleton at rates dictated by the local chemistry. As the two ends 
of a cytoskeletal filament loaded with motors carrying cargo are potentially 
well separated in space, their chemical environments need not be identical. It 
is thus possible that motor loading and unloading rates could be different at 
either end of the filament. The generalization of ratchet models for the motion 
of interacting motors, extended to allow for open boundary conditions, could 
potentially allow for the boundary injection and removal of motors. This in- 
jection and removal will compete, in general, with the bulk equilibration via 
Langmuir kinetics[34,35,36] of the motor density along the filament, but we 
will ignore the possibility that motors are lost or gained along the filament, 
accounting only for their entry and exit at the boundaries. For recent work 
which incorporates both bulk non-conservation and the effects of geometry, 
see Refs. [37,38,39,40,41]. 

For simulations with open boundary conditions, we attach two boundary sites 
to the open chain with A^ sites; the total number of sites is then A^ + 2. These 
boundary sites are filled with probability a (at the right boundary) or /?(at the 
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r = 0.05 








Fig. 8. Current versus the boundary and injection rates a and (3, obtained through 
Monte Carlo simulations for parameter values N = Lw = 128 * 10 sites, r = 0.05, 
Pio = 0.04 and Pqi = 0.05. 

left boundary). Hopping rules at these special sites are chosen to be consistent 
with the direction of the net current in the absence of boundary driving and 
motors at these singled-out sites always hop unidirectionally. Results for this 
model are shown in Fig 8 which exhibits the current in the system J as a 
function of the input and output rates, a and 13. Note the striking feature 
of the plot - the remarkable independence of J on a and 13 for a wide range 
of these parameters. For this range of input and output rates, the current is 
not only independent of the rates at which motors are added or subtracted at 
both ends but the current being passed through the system is pegged at its 
maximum value. Fig. 9 shows the steady state densities at varying input and 
output rates. 

Note the existence of three distinct phases: (I) a phase in which the current 
and density can be varied by changing only j3 and is independent of a, (II) a 
phase in which the current and density can be varied by changing only a and is 
independent of /3, (III) a "constant current" phase in which the current attains 
its maximum possible value and is independent of both a and f3. The boundary 
driven ASEP has a phase diagram which is qualitatively very closely similar, 
with a low density phase, a high density phase and a intermediate constant 
current phase[20]. It should, in principle be possible to investigate the same 
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Fig. 9. Current versus the boundary and injection rates a and f3, obtained through 
Monte Carlo simulations for parameter values N = Lw = 128 * 10 sites, r = 0.05, 
Pio = 0.04 and Pqi = 0.05. 

effects described in Refs. [34,35,36], where Langmuir kinetics competes with 
one-dimensional transport to yield a phase diagram with complex structure. 
We do not, however, address this interesting problem here. 



4 Motor Microtubule Pattern Formation 



The mitotic spindle of a dividing eukaryotic cell is a self-organized struc- 
ture at the sub-cellular scale. Such structures are usefully thought of as pat- 
terns, by which we mean spatially inhomogeneous yet stable steady states, 
defined through the interaction of motors and microtubules. Experiments on 
centrosome-free fragments of the cytosol containing both motors and micro- 
tubules obtain self-organized radial structures called asters. Single asters, in 
addition to other complex patterns such as vortices, disordered aster-vortex 
mixtures and lattices of asters and vortices, are also seen in vitro, in exper- 
iments on mixtures of molecular motors and microtubules in a quasi-two- 
dimensional geometry [21]. 

What physical processes stabilize such structures? Direct molecular dynamics 
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simulations which incorporate any level of molecular-scale realism are cur- 
rently incapable of tacking the pattern formation problem. One must then 
rely on approximate models for such systems, again keeping in mind the ne- 
cessity of maintaining contact with biological reality: our models must be as 
simple as possible but no simpler. This problem has been studied extensively 
over the past 5 years or so with significant contributions from several groups 
[44,48,45,46,47]. This section summarizes work in this direction first presented 
in Ref. [7]. 

In this section, hydrodynamic equations of motion for a coarse-grained field 
representing the local orientation of microtubules as well as for local motor 
density fields are described. Our model treats motors attached to microtubules 
differently from motors which diffuse freely in solution. Motors which move 
on microtubules are referred to as "bound" motors, while those which diffuse 
in the ambient solvent are referred to as "free" motors. These are described 
by coarse-grained fields denoted by rrib and rrif respectively and obey different 
equations of motion. 

We will take microtubules to be oriented by complexes of bound motors, yield- 
ing patterns at large scales. Our results are: We obtain a single vortex as a 
stable final state for large motor densities in some regimes. However, in other 
regimes, asters are favoured. A "lattice of asters" state is stabililized in our 
model through a low-order relevant term in the equation of motion for the 
microtubule orientation. On small systems, constraints due to confinement 
favour a small number of asters, whose number can be increased systemat- 
ically as parameters are varied. We have also calculated the distribution of 
free and bound motors in asters and vortices obtained in our model; we derive 
an exponential decay of bound motor densities away from aster cores, mod- 
ulated by a power-law in which the exponent of the power law depends in a 
non-universal way on dynamical parameters. 

In the absence of inter conversion terms changing a bound motor to a free 
motor, mb obeys a continuity equation involving the current of motors trans- 
ported along the microtubules. The free motor field mj obeys a diffusion 
equation with a diffusion constant D. These two fields are coupled through 
mechanisms which convert "free" motors to "bound" motors and vice versa. 
We will take 7j^b and 75_+ / to be the rates at which free motors become bound 
motors ("on" rate) and vice- versa ('off" rate). 

Our aim is to write down an minimal set of equations capable of both de- 
scribing the variety of patterns formed in such interacting motor-microtubule 
mixtures. We are guided both by symmetry, as is appropriate for a fully non- 
equilibrium system in which detailed balance based on rates derivable from a 
hamiltonian does not exist, as well as by considerations of simplicity: of an 
infinity of possible non-equilibrium terms allowed in our equations of motion. 
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we choose the simplest. In appropriately scaled units, the equations then are 
dtruf = V^nif - ^f^b^rif + ^b^jmh (8) 

dtrrib = -V ■ (m^T) + ^f^trrif - %^fmb (9) 

dtT = CT(1 - T^) + rribV^T + eVm^ ■ VT + /tV^T + SVrrib (10) 

Free and bound motor density profiles in vortex and aster configurations can 
be obtained from the above equations. We set the time derivatives to zero i.e. 
dtiTLf = dtiTLh = 0, obtaining 

V^ruf - V ■ (mfoT) = 0, (11) 

V^m/ + {'yb-^firib - 'jf-.b'mf) = 0. (12) 

For a single vortex i.e. T = 6, we may assume radial symmetry and thus 
rrif = Tnf{r) and rUb = mb{r). We then obtain, 

m/(r) = ci + C2ln(r), (13) 

where ci and C2 are constants to be determined by boundary conditions and 
normalization. The relation '~fb-*fiTib — 7/_»bm/ = yields 

^6(r) = ^^^^(ci + C2ln(r)). (14) 

For the motor distribution about a single aster i.e. T = — f, we again as- 
sume radial symmetry for the bound and free motor densities. The boundary 
condition that the total motor current vanishes at the boundary implies 

drmf{r) = —mb{r). (15) 

Thus we have 

I 
dlnif + ( 'jb^f)drmf - 'jf^birif = 0. (16) 

The general solution to the equation above is a combination of confiuent hy- 
pergeometric functions and can be written in terms of the two solutions of the 
hypergeometric Kummer equation, modulated by an exponential. It is useful 
to define a quantity p given by 

P = lil- i,^'^\ )■ (17) 

2 y^7L/ + 47/^6 



Note that < p < 0.5 with 7^^/, 7/^5 > 0. 

The asymptotics is obtained using an integral representation for the appro- 
priate hypergeometric functions. Our final result for the motor distribution 
about a fixed aster configuration is 



.-rli 






(7L/ + ^V^i^r r-P 



(7Lf + 47/^.r/VV e^' 



with r^ = l ^^'^^ V(^^^/+ i2Z^| = |lg_^|. The correlation length ^ and the 
power-law exponent p depend on 7/-+6 and 76^/. We see that the bound motor 
density in the aster case has an exponential fall modulated by a power-law 
tail. 

We also solve Equations 8,9 and 10 numerically on an L x L square grid indexed 
by {i, j) with i = 1, . . . L and j = 1, . . . L. The equations for the free and bound 
motor densities are evolved using an Euler scheme. We impose the boundary 
condition that no current (either of free or bound motors) fiows into or out of 
the system. This condition is easily imposed by setting the appropriate current 
to zero. The T equation is differenced through the Alternate Direction Implicit 
(ADI) operator splitting method in the Crank-Nicholson scheme. 

Our simulations are on lattices of several sizes, ranging from L = 30 to 
L = 200. We vary the motor density in the range 0.01 to 5 in appropriate 
dimensionless units. We work with two different types of boundary conditions 
on the T field. In the first, which we refer to as refiecting boundary conditions, 
the microtubule configuration at the boundary sites is fixed to point along the 
inward normal. In the second, which we refer to as parallel boundary con- 
ditions, microtubule orientations at the boundary are taken to be tangential 
to the boundary. In both these sets of boundary conditions, the state of the 
boundary T vectors is fixed and does not evolve. The total number of mo- 
tors, initially divided equally between free and bound states and distributed 
randomly among the sites, is explicitly conserved. 

Figures 10(a) - (d) depict four stable configurations obtained in different 
regimes of parameter space for an L = 100 lattice. Fig. 10(a) shows a dis- 
ordered arrangement of microtubules obtained at very low motor densities 
{m = 0.005) with e = 0.5 and S = 0. Figure 10(b) shows an aster-vortex 
mixture obtained at m = 0.01 at the same values of e and S. This figure 
is to be contrasted to Fig. 10(c), obtained at m = 0.05, taking e = 5 and 
S = 0.001. Note the absence of asters in this regime of parameter space. Fi- 
nally, Fig 10(d), obtained with m = 0.5, e = 1 and 5* = 1, illustrates a lattice 
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Fig. 10. Steady state configurations in our model at different parameter values (see 
text): (a) Disordered states obtained at very low motor densities [m = 0.005 e = 0.5 
and S" = 0]; (b) Aster-vortex mixture obtained at [m = 0.01, e = 0.5 and S* = 0]; 
(c) Lattice of vortices at [m = 0.05, e = 5 and S = 0.001]; (d) Lattice of asters 
obtained at [m = 0.5, e = 1 and S* = 1] 

of asters, with asters being the only stable defects present. We can vary the 
sizes and numbers of asters obtained in configurations such as the one shown 
in Fig. 10(d), by changing S. A larger 5* yields a large number of small asters, 
while smaller values of S yield a smaller number of large asters [24]. 

Our results are summarized in Fig. 11 which shows the states which dominate 
in the three-dimensional space spanned by e, S and m. For S* = 0, we obtain 
disordered/aster-vortex mixture states at low motor density, which become a 
lattice of vortices at somewhat higher motor densities. Large values of e (e > 
1) yield well-formed vortex-like configurations while small e yields structures 
better described as aster-vortex mixtures. At intermediate values of e and m, 
spirals rather than vortices appear to dominate. At large m, with S = and 
large e, a single vortex is obtained[44,48]. 
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Fig. 11. Qualitative map of steady states illustrating how different states, the 
disordered state, the aster-vortex mixture state, the lattice of vortices state, the 
single vortex, the lattice of asters, dominate in different regimes of parameter space; 
for a definition of parameters see text. The parameter e is plotted on the y axis, with 
the total motor density m plotted on the x axis. The parameter S extends out of 
the e — ra plane. Of the states shown, the lattice of asters is obtained generically for 
non-zero S (out of the plane of the figure) , whereas the other states are associated 
with the S" = plane, although they appear to survive provided S is small enough. 

For non-zero but small S*, these states appear to continue out of the S* = 
plane but are rapidly replaced by a lattice of asters for larger S. A cut of 
Fig. 11 at finite S yields disordered states at small m and a lattice of asters 
at larger m. We can thus understand the sequence of patterns formed upon 
increasing m in mixtures of kinesin constructs with microtubules in terms of 
a trajectory which begins in the S* = (or S sufficiently small) plane in the 
disordered phase and transits between the aster-vortex mixture and the lattice 
of vortices (both of which lie in this plane) as m is increased. As m increases 
further and the effects of the S term become important, such a trajectory 
moves out towards non-zero 5, encountering the lattice of asters. 

We have also examined the effects of changing the motor processivity, a quan- 
tity proportional to the ratio of 7/^6 to 76^/. Smaller values of this ratio are 
appropriate to molecular motors such as NCD. At 7/-^^ = 0.005, 7b_^/ = 0.05, 
we find that the disordered regime shown in Fig. 1 expands, so that at equiv- 
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alent values of m disordered states occupy much of the domain associated 
previously with the lattice of vortices. Whereas kinesins follow the sequence 
disordered - lattice of vortices - aster vortex mixture - lattice of asters as 
the m is increased, a mixture of microtubules with NCD motors bypasses the 
lattice of vortices altogether, transiting directly from the disordered state to 
the lattice of asters in the experiments [22]. In terms of Fig. 11, the expanded 
regime of disordered states for the NCD motor suggests that patterns such as 
the lattice of vortices and the aster- vortex mixture may be inaccessible at the 
motor densities at which the experiments are done, since the effects of the S 
term might be expected to dominate at large m. 

It is interesting to note that the generic state which is obtained at large m is 
a lattice of asters for non-zero S, in contrast to the predictions of the earlier 
model of Lee and Kardar, which indicates that the large motor density state 
is always a single vortex. This feature, a direct consequence of the presence 
of the crucial SVm term, agrees with experiment. Further results, including 
a discussion of motor-motor interactions in the context of motor- microtubule 
pattern formation, qualitative "free-energy" based arguments for the stability 
of patterns and a detailed discussion of the effects of confinement can be found 
in Ref. [7]. 



5 Conclusions 



The interacting thermal ratchet model discussed in the first part of this paper 
introduces a primitive level of biological realism into models for the motion 
of interacting motor proteins. While the model itself may not appear appre- 
ciably more realistic than the ASEP itself, the use of ratchet models gener- 
alized to include interactions is attractive, since such models incorporate the 
true reason for the symmetry breaking which occurs when molecular motors 
move unidirectionally in a Brownian environment. This reasoning is obscured 
in the currently popular ASEP-based models, where the symmetry breaking 
which leads to a non-trivial current is imposed by hand, through the defi- 
nition of the hopping rates. Making the connection to ratchet models also 
facilitates the understanding of many issues relating to the coherence or re- 
liability of transport [49] in motor systems as well as of efficiency in energy 
transduction [50]. The exclusion processes simply lack the necessary structure 
for such discussions to be meaningful. 

We have also described calculations which relate to a non-equilibrium pattern 
formation problem: the formation of patterns in mixtures of microtubules and 
molecular motor constructs. The model is able to reproduce virtually all the 
patterns obtained in the experiments, but has the advantage over direct simu- 
lations that the number of parameters which need to be included is small. We 
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have also been able to derive many of the features of the patterns which form, 
including the sequence of patterns which are obtained as the motor density is 
increased. Our model rationalizes several features of the experiments, in many 
cases for the first time. These include: (i) the sequence of patterns obtained as 
the motor density is increased, (ii) the prevalence of the lattice of asters state 
and (iii) the difference in the sequence of patterns formed by conventional 
kinesins and the NCD motor. Many further features of these equations are 
currently being explored. 



One final point relates to the nature of steady states which are obtained in the 
ratchet models generalized to include interactions between motors and has to 
do with the possibility of "robustness" in these models. Barkai and Leibler[42] 
and Alon et. al. [43] study the chemotaxis network of E. Coli, suggesting that 
this network exhibits exact and robust adaptation, over a wide range of vari- 
ation of parameters. The robustness of adaptation in this case is an example 
of biochemical robustness, since it has its origins in specific features of the 
biochemical network underlying chemotaxis in E. Coli. The "tensegrity" of 
the cytosketal network of living cells can be thought of as another form of ro- 
bustness. Such robustness of the structural elements in the cell to mechanical 
perturbations is distinct from biochemical robustness. It is thus interesting 
to ask the following question: Are other manifestations of robustness possible 
and are there biological situations in which they may be relevant? 



We suggest another intriguing possibility for robust behaviour in biological 
systems, the robustness of certain non- equilibrium steady states of biological 
systems to wide classes of physical perturbations. A simple illustrative exam- 
ple in the context of the boundary driven interacting ratchet model is the 
independence of the steady state current on a and j3 for a very large range of 
such parameter values. This indicates that the current is insensitive to fairly 
large fluctuations in the input and output rates and is, moreoever, pegged 
to the largest possible value it can take. The robustness here is for physical 
reasons, essentially having to do with the one-dimensional character of this 
steady state and the fact that the system can adjust its density profile to ac- 
comodate changes in the boundary driving rates. While allowing for Langmuir 
kinetics in the bulk in this specific case will generically destroy such behaviour 
in the thermodynamic limit, small systems should exhibit the same qualitative 
behaviour. It remains to be seen if such "steady-state robustness" or "physical 
robustness" is in fact a feature of some aspects of in vivo cellular function, 
independent of the models we contrive to describe such function. It would be 
interesting to look at other examples of non-equilibrium steady states in bio- 
logical systems, possibly those which involve cell-scale flows (as, for example, 
in cell streaming), to see if this idea might find support. 



23 



6 Acknowledgements 



The work described in the first section was done in collaboration with Yashar 
Aghababaie and Michael Plishke at Simon Fraser University, with research 
supported by the NSERC of Canada. The work described in the third part, 
relating to the motor-microtubule pattern formation problem was done in 
collaboration with Sumithra Sankararaman at the Institute of Mathematical 
Sciences and P.B. Sunil Kumar at the IIT, Madras and was partially supported 
by the DST (India). I am grateful to these collaborators for all I have learnt 
from them. At various stages, I have benefited from conversations with Michael 
Wortis, Jacques Prost, Mustansir Barma, Deepak Dhar, Arun Jayannavar, 
Madan Rao, Sriram Ramaswamy and Amit Kumar Bhattacharjee. I would 
also like to thank Debashish Chowdhury for organizing a very stimulating and 
topical conference on these and related issues where this work was presented. 



References 



[I] S. Ramaswamy and M. Rao, C.R. Acad. Sci Ser IV, Phys. Astrophys. 2, 818 
(2001) 

[2] Y. Hatwalne, S. Ramaswamy, M. Rao and R. A. Simha, Phys. Rev. Lett., 92 
118101 (2004) 

[3] J. Toner, Y. Tu and S. Ramaswamy, Annals of Physics, 318, 170 (2005) 

[4] A.J. Koch and H. Meinhardt, Reviews of Modern Physics, 66, 1481 (1994). 

[5] D. Chowdhury, A. Schadschneider and K. Nishinari, Physics of Life Reviews, 
2, 318 (2005). 

[6] Y. Aghababaie, G.I. Menon and M. Plischke, Phys. Rev. E, 59, 2578 (1999). 

[7] S. Sankararaman, G. I. Menon and P. B. Sunil Kumar, Phys. Rev. E, 70, 
031905 (2004) 

[8] A.D. Mehta, M. Rief, J. A. Spudich, D.A. Smith and R.M. Simmons, Science, 
283, 1689 (1999) 

[9] J. Howard, Mechanics of Motor Proteins and the Cytoskeleton, Sinauer, 
(Sunderland, MA, U.S.A, 2001). 

[10] F. Jiilicher, A. Ajdari and J. Prost, Rev. Mod. Phys. 69, 1269 (1997). 

[II] P. Reimann, Phys. Rep., 361, 57, (2002). 

[12] F. Jiilicher and J. Prost, Phys. Rev. Lett, 75, 2618 (1995) 
[13] F. Jiilicher and J. Prost, Phys. Rev. Lett, 78, 4510 (1997) 



24 



[14] S. Camalet, T. Duke, F. Jiilicher and J., Prost, Proc. Nat. Acad. Sci USA, 97 
3183 (2000) 

[15] D. Dan, A. M. Jayannavar and G. I. Menon Physica A, 318, 40 (2003) 

[16] I. Derenyi and T. Vicsek, Phys. Rev. Lett. 75, 374 (1995); I. Derenyi and A. 
Ajdari, Phys. Rev. E 54, R5 (1996). 

[17] K. Nishinari, Y. Okada, A. Schadschneider and D. Chowdhury, Phys. Rev. 
Lett, 95 118101 (2005) 

[18] J. Krug, Phys. Rev. Lett, 67 1882 (1991) 

[19] B. Derrida, E. Domany and D. Mukamel, J. Stat. Phys, 69 667 (1992); G. 
Schutz and E. Domany, J. Stat. Phys, 72 277 (1993); B. Derrida et. al. J. 
Phys. A, 26 1493 (1993). 

[20] R.B. Stinchcombe, Advances in Physics 50 431 (2001). 

[21] F. Nedelec, T. Surrey, A. C. Maggs and S. Leibler, Nature 389, 305 (1997). 

[22] T. Surrey, F. Nedelec, S. Leibler and E. Karsenti, Science 292, 1167 (2001). 

[23] F. Nedelec and T. Surrey, C. R. Acad. Sci. Paris, Serie IV, 841 (2001). 

[24] F. Nedelec, T. Surrey, A. Maggs, Phys. Rev. Lett 86, 3192 (2001). 

[25] F. Nedelec, .Journal of Cell Biology, 158, 1005 (2002). 

[26] F. Nedelec, T. Surrey and E. Karsenti, Current Opinion in Cell Biology 23, 
118 (2003). 

[27] R. Heald, R. Tournebize, A. Habermann, E. Karsenti, and T. Hyman, 
Jour. Cell Biol. 138, 615 (1997); R. Heald, R. Tournebize, T. Blank, R. 
Sandaltzopoulos, A. Hyman, and E. Karsenti, Nature 382, 420 (1996). 

[28] M. Kardar, G. Parisi and Y. C. Zhang, Phys. Rev. Lett. 56, 889 (1986). 

[29] See A.-L. Barabasi and H.E. Stanley, Fractal Concepts in Surface Crowth, 
(Cambridge University Press, 1995) for a review of much of this work, 
especially in the context of nonequilibrium surface growth. 

[30] H. van Beijeren, R. Kutner and H. Spohn, Phys. Rev. Lett. 54, 2026 (1985). 

[31] D. Dhar, Phase Transitions, 9, 51 (1987). 

[32] L-H. Gwa and H. Spohn, Phys. Rev. A, 46, 844 (1992). 

[33] M. Plischke, Z Racz and D. Liu, Phys. Rev. B 35, 3485 (1987). 

[34] A. Parmeggiani, T. Franosch and E. Frey, Phys. Rev. Lett, 90, 86601 (2003) 

[35] A. Parmeggiani, T. Franosch and E. Frey, Phys. Rev. E, 70, 46101 (2004) 

[36] H. Hinsch, R. Kouyos and E. Frey, cond-mat/051244 '^ 



25 



[37] R. Lipowsky, S. Klumpp and T.M. Nieuwenhuizen, Phys. Rev. Lett, 87, 108101, 
(2001) 

[38] T.M. Nieuwenhuizen, S. Klumpp and R. Lipowsky, Europhys. Lett, 58, 468 
(2002). 

[39] T.M. Nieuwenhuizen, S. Klumpp and R. Lipowsky, Phys. Rev. E, 69, 061911 
(2004). 

[40] S. Klumpp, A. Mielke and C. Wald, Phys Rev E, 63, 031914 (2001). 

[41] S. Klumpp, T.M. Nieuwenhuizen and R. Lipowsky, Physica E, 29, 380 (2005). 

[42] N. Barkai and S Leibler, Nature (London), 387, 913 (1997) 

[43] U. Alon, M.G. Surette, N. Barkai and S. Leibler, Nature, (London), 397 168 
(1998) 

[44] H. Y. Lee and M. Kardar, Phys. Rev. E 64, 056113 (2001). 

[45] T. B. Liverpool and M. C. Marchetti, Phys. Rev. Lett, 90, 138102 (2003); F. 
Ziebert and W. Zimmermann, ibid. 93, 159801 (2004). 

[46] K. Kruse et al, Phys. Rev. Lett 92, 078101 (2004) 

[47] I. S. Aranson and L. S. Tsimring, Phys. Rev. E, 71, 050901(R) (2005) 

[48] J. Kim, Y. Park, B. Kahng, and H.Y. Lee, Jour. Kor. Phys. Soc. 42, 162 
(2003). 

[49] J. A. Freund and L. Schimansky-Geier, Phys. Rev. E, 60 1304 (1999); T. Harms 
and R. Lipowsky, Phys. Rev. Lett, 79, 2895 (1997); R. Krishnan, D. Dan and 
A.M. Jayannavar, Physica A, 354, 171 (2005) 

[50] A. Parmeggiani, F. Jiilicher, A. Ajdari and J. Prost, Phys Rev. E, 60, 2127 
(1999) 



26 



